
## look at change in graduation probability
grad<-read.csv("./estimation_stuff/grad.csv")
tmp<-merge(student_data, grad, by.x="pid", by.y="id")
tmp$avg_gpa<-with(tmp, (gpa1 + gpa2)/2)

y_model_baseline_avg_df<-ddply(results_baseline_both_sem[,c("pid", "y_model")], .(pid), summarise, avg_y_model_baseline=mean(y_model, na.rm=TRUE))
y_model_cf_avg_df<-ddply(results_cf_both_sem[,c("pid", "y_model")], .(pid), summarise, avg_y_model_cf=mean(y_model, na.rm=TRUE))
y_diff_avg_df<-ddply(results_diff_both_sem, .(pid), summarise, avg_y_diff=mean(y_diff, na.rm=TRUE))
tmp<-merge(tmp, y_model_baseline_avg_df, by="pid")
tmp<-merge(tmp, y_model_cf_avg_df, by="pid")
tmp<-merge(tmp, y_diff_avg_df, by="pid")
tmp<-merge(tmp, subset(results_baseline_long_sem1, select=c(pid, mu_y)), by="pid")

grad_probit<-glm(grad ~ mu_y + avg_y_model_baseline, family=binomial(link="probit"), data=tmp, na.action="na.exclude")
grad_probit_vars<-glm(grad ~ black + male + hsgpa +  combact + studyhs + estudy  + avg_y_model_baseline, family=binomial(link="probit"), data=tmp, na.action="na.exclude")

pred_grad_baseline<-predict(grad_probit, type="response")
pred_grad_cf<-predict(grad_probit, newdata = data.frame(mu_y=tmp$mu_y,avg_y_model_baseline = tmp$avg_y_model_cf),type="response")
diff_pred<-pred_grad_cf - pred_grad_baseline

pred_grad_df<-data.frame("base_grad_prob" = t(cbind(total=mean(tmp$grad, na.rm=T), t(tapply(tmp$grad, tmp$race, mean, na.rm=T)), t(tapply(tmp$grad, tmp$sex, mean, na.rm=T)), t(tapply(tmp$grad, tmp$hsgpa_level, mean, na.rm=T)))), 
                         "change_pred_grad" = t(cbind(total=mean(diff_pred, na.rm=T), t(tapply(diff_pred, tmp$race, mean, na.rm=T)), t(tapply(diff_pred, tmp$sex, mean, na.rm=T)), t(tapply(diff_pred, tmp$hsgpa_level, mean, na.rm=T)))))  
pred_grad_df$cf_grad_prob<-pred_grad_df$base_grad_prob + pred_grad_df$change_pred_grad

pred_grad_df<-pred_grad_df[,c(1,3,2)]
print(xtable(pred_grad_df, caption="Change in predicted graduation rate", label="tab:change_pred_grad", digits=3), file="./output/change_pred_grad.tex", caption.placement ="top")

stargazer(grad_probit, grad_probit_vars, out="./output/grad_probit.tex", title="Probit of graduating on average of first year GPA", label="tab:grad_probit")

